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The dynamics of the triple gas-liquid-solid contact line is analysed for the case where 
the gas is the saturated vapour corresponding to the liquid, like in the vapour bubble in 
boiling. It is shown that even small superheating (with respect to the saturation tem- 
perature) causes evaporation of the adsorption liquid film and the true triple contact 
is established. It is shown that the hydrodynamic contact line singularity cannot be 
relaxed with the Navier slip condition under such circumstances. Augmented with the 
second derivative slip condition is proposed to be applied. For the partial wetting condi- 
tions, a non-stationary contact line problem where the contact line motion is caused by 
evaporation or condensation is treated in the lubrication approximation in the vicinity 
of the contact line. High heat fluxes in this region require the transient heat conduction 
inside the heater to be accounted for. Two 2D problems, those of drop retraction with 
no phase change and of drop evaporation are solved and analysed as illustrations of the 
proposed approach. 
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1. Introduction 



The modeling of bubble growth remains to be poorly understood in spite of its impor- 
tance for estimation of the overall heat transfer during the nucleate boiling. It is widely 
recognised that a large part of the boiling heat transfer is due to the evaporation at the 
foot of the bubble, in the region often called "microregion" . Therefore a poor knowledge 
of the microregion processes can lead to a large uncertainty of the overall heat transfer 
estimations. 

One of the most important advanc es in the understanding of this problem and of the 



boiling in general has been made by 



Cooper fc Llovdl (|l969| ) when they suggested the 



existence of a thin liquid layer ( "microlayer" ) that separates the bubble from the heater. 
They have based their conclusion on the hcmisphcricity of the bubbles. As a matter of 
fact, the bubbles become hemispherical under the action of the inertial hydrodynamic 
forces which are important when the bubbles grow fast enough, that is at low system 
pressureflj. The present understanding of the study of the micr olayer evaporation is based 



on the approach developed originally bv lWavner et al 



(|1976f ) for the evaporation of the 
continuous liquid meniscus. Among numerous a pplications of these ideas in boiling and 
evaporation in general one can cite the works (jStephan fc Hammer 



1999; 


Morris 


2000; 


Aiaev et al. 


2002) 



1985 



Son fc Dhir 



assumed that the thinnest part of the microlayer ( "adsorption film" ) does not evaporate 
due to the van der Waals attraction forces that exist between the molecules of the solid 
heater and the fluid. 

While the theoretical and numerical studies deal with the continuous microlayer, ex- 
f At high pressures (i.e. comparable with the critical pressure for the given fluid), the bubble 

shape is rather spherical than hemispherical. Cooper and Lloyd were cautious about applicability 

of their idea to this case. 
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van Ouwerkerk 



Nikolavev et al. 



1972; 



Torikai et al. 


1991 


TheoDhanous et al. 


2002 



Hegs eth et al 



2005 



2006) when the heat flux transferred from the heater to the fluid is com- 
parable to the critical heat flux (CHF) at which the heater dries out and the transition 
to the film boiling takes place. This phenomenon is called "boiling crisis" or "departure 
from nuclea te boiling" . To predict the CHF, the understand ing of the drying dynamics 



is essential (jNikolavev fe; Bevsens 



1999; 



Nikolavev et al 



2001). It is likely that dry spots 



exist at much smaller fluxes however they are smaller and thus more difficult to observe. 

The drying implies the direct contact of the vapour with the solid and the existence 
of the triple vapour-liquid-solid contact line (CL). The determination of the CL position 
becomes a key to the dryout simulation and requires the clear understanding of the heat 
transfer near the CL, which is the subject of this article. In the next section we discuss 
the conditions under which the CL concept is relevant. 

Besides the description of drying, the CL issue is important for the understanding of 
the bubble departure from the heater, another key feature of boiling. The surface tension 
is the only force that provid es the adhesion of the bubble at the last, thermally controlled 



stage of the bubble growth (jSon fc Dhir 



19991 ). This force exists only when the CL exists 



and is proportional to its length. The CL position is thus needed to be known to assess 
the bubble departure size. 

High heat fluxes that occur in the microregion show the importance of the accounting 
for the thermal cond uction in the heater tha t influences strongly the evaporation rate at 



the bubble interface (jNikolavev et al. 



2001). It has not been taken into account in the 



previous studies of the microregion and will be discussed in the present work. 

The hydrodynamic stresses are well known to be infinite at the moving CL when the 
classical no-slip boundary condition (zero tangential velocity at the solid wall) is applied. 
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This non-physical singularity need to be relaxed. The replacement of the no-slip condition 
by the Navier slip condition presents a singularity relaxation solution that is used by most 
researchers working in this field. In sec. 14.11 we show that this solution in inappropriate 
when heat transfer is taken into account. We propose a new slip boundary condition that 
allows the CL singularity to be relaxed self-consistently. 

Two coupled parts of the CL problem can be distinguished: the heat transfer problem 
which is dealt with in the section [3] and the determination of the shape of the liquid- 
vapour interface near the CL discussed in the section [4j However one needs to analyse 
first the conditions at which the liquid film at the heater dries out and the CL actually 
appears. 



2. Does the adsorption liquid film exist at the heater? 

Numerous experimental studies have been carried out to detect the dry spots on the 
hot enough heater. In the experiments 



van Ouwerkerk 



1972 



Torikai et al. 



1991D with a 



transparent heater the total internal reflection at the heater-fluid interface was used to 
detect the dry spots. This shows that even if the liquid film existed, its thickness was 
smaller than the about one tenth of the light wavelength (~ 50 nm) . In the experimenta l 
study of evaporation at extremely high (nearly critical) pressures (|Hegseth et aLl 120051 ) . 
the motion of the CL separating the wetted and dry areas have been studied optically, 
with about the same precision. In order to check if the thinner fil m exist, only molecular 



dynamics simulations data is available. Most simulation studies ([Maruvama fc Kimura 



2000: 



Consolini et al 



20031 ) show only separate molecules adhering to the solid so that the 



co ntinuous liqui d layer is nonexistent. The liquid layer was rather thick in the simulations 



of 



Freundl (|2005l ) but only the case of almost insignificant superheating have been studied. 



In the absence of the other data, let us analyse this problem theoretically. The studies 
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of thin films evaporation analogous to those by Wayner belong to the class of the "mean 

field" theories. This term comes from the statistical physics and means that the thermal 

fluctuations are neglected. It is well known that the fluctuations can be neglected for large 

systems. In the present context, "large" means that the number of the molecular layers 

in the film is large. Let us estimate the film thickness for which the thermal fluctuations 

need to be taken into account. 

The van der Waals energy of a fluid molecule situ ated at a distance d from the half- 



space filled with the material j can be calculated as ([Israelachvili 



1992) 



W {d) = -B Fo n 3 — , (2.1) 

where nj means the number of molecules per unit volume, S or L should be used instead 
of j depending on if the half-space is filled by solid (i.e. heater) or liquid respectively, B 
is the constant of r~ 6 attraction between respective molecules. Note that the expression 
for the disjoining pressure pa = A/d 3 can be calculated from (|2.ip . 

It is evident that the van der Waals forces cannot hold the fluid molecules at any 
superheating: at some temperature the kinetic energy of the molecules (usually neglected 
in the previous analyses of the adsorption films) will become so high that they will leave 
the heater. In order to find at which superheating the adsorbed liquid film will disappear, 
one needs to compare the cohesion and thermal (kinetic) energies for a molecule sitting 
at the vapour-liquid interface. Compare two configurations (Fig. [!}. The first is the 
plane vapour-liquid interface at saturation temperature T sat (Fig. [TJi). In the second 
configuration, a part of the lower half-plane situated at the distance d from the interface 
is replaced by a solid heater so that a liquid film is formed (Fig. [TJd). The change of 
cohesion energy of an interface probe molecule between two these configurations is 

AW(d) = W L (d) - W s {d) = {B FS n s - B FF n L )-^ (2.2) 
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(a) (b) 
Figure 1. To the adsorption film estimation. 
The excess kinetic energy gained by this probe molecule due to the larger (for the second 
configuration) interface temperature T % is ksiT 1 — T sat ), where ks is the Boltzmann 
constant. By solving the equation AW(d) = fcs(T 4 — T sat ) for d, one obtains the critical 
thickness d c . The molecules situated at larger distances are not kept by the van der Waals 
forces in the film and are evaporated. 

To obtain the upper bound for d c , we assume the extreme case Bps n s 3> BffTIl 
that strongly favou rs the wetting film formation. By using the order of value estimations 



( Israelachvili 



1992h n s ~ 10 29 m~ 3 and Bps ~ 1CT 77 J-m 6 , one obtains 4 = 3 nm for the 
superheating T l — T sat = 1°C. Since for such thin films the temperatures of the interface 
and of the heater almost coincide and the local superheating attains often muc h larger 



value s, d c value in reality is even smaller. Such films adhere rather than flow (jFreimr 



20051) . As far as the fluid flow is concerned, their existence can thus be neglected. The 
conclusion is that even small superheating causes complete evaporation of the liquid film 
and formation of the CL. Therefore the classical theory of meniscus evaporation from 
which the CL is absent (the "interline" was defined in these models as a contact between 
the vapour, bulk liquid and the film) needs to be revisited. The term "microregion" will 



be used below to denote the vicinity of the CL. 




solid heater 



Figure 2. Geometry of the general problem. The chosen direction of the normal to the 

interface is shown. 

When the CL is present, one cannot use the diverging at d — ► expression for pd- 
It must be corrected for the short range intermolecular interac tions to describe self- 



consistently the partial wetting case (jBrochard-Wvart et al. 



1991J). However the disjoin- 



ing pressure ceases to be an important factor for the wedge geometry. Indeed, the area 
where pd is significant (h < 100 nm) is very small. Followin g previous works on th e 



treatment of heat transfer in the case where the CL is present ((Anderson &: Davis 
we thus omit pd from further consideration. 



19951) . 



3. Heat transfer in the microregion 

In the description of the microregion, the curvature of the vapour-liquid interface in 
the direction parallel to the CL can be neglected with respect to the curvature K in the 
perpendicular direction. The interface can thus be described by its 2D contour (Fig. [2]) . 
We shall assume that the contact angle 9 ^ is small (which is generally the case for most 
heaters) so that the slope of the interface is small. The interface can then be described by 
its hight y = h(x) and the small slope assumption means \dh/dx\ <C 1. All the variables 
can then be described as functions of x. 
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3.1. Problem statement in a general case 

The most naive approach to the contact line heat transfer would be assuming both 
the liquid-vapour and liquid-heater interfaces isothermal. The background argument is 
that the temperature of the vapour-liquid interface T l is generally very homogeneous 



( Marek fe Straub 



2001 



Barthes et al 



20071 ) (T z = T sat ) in the absence of impurities 



in both liquid and vapour phases (which is the case considered here). On the other 
hand, the surface of the metal heater is often assumed isothermal due to its high thermal 
conductivity. To vaporise the liquid, the heater surface temperature T$ needs to be higher 
than T sat . Since the CL belongs to the both interfaces, its temperature is ambiguous, 
which is well known to generate the non-integrable divergence of the heat flux at the 
vapour liquid interface (assumed to be positive at evaporation) q l L (x) ~ 1/x, which 
means that the integral heat balance cannot be satisfied at all. To avoid this difficulty, 
the temperature needs to be allowed to vary along at least one of the interfaces. 

T % can only vary due to kinetic effects that are tak en into account by introduction of 



the interface resistance R z ([Stephan fe Hammerl ll985) 



where 



R l = 



T i = T? at + R i qi, 



T sat y/2TrR g T sa t/ M(pl - pv) 
H 2 p L pv 



(3.1) 



M is the molar weight, R g is the universal gas constant, H is the latent heat and pl 
(pv) is the liquid (vapour) density. Rather than the saturation temperature for the given 
system pressure T sat , (|3.1[) involves the saturation temperature 

Ap 



rpp _ rp 

± sat — A sat 



that depends (jWavner et al 



1976; 



Mm. 



(3.2) 



Stephan fc Hammei 



1985) on the pressure jump Ap 



Pv^Pl across the interface. Ap includes contributions of all forces that act on the system. 
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If other forces like gravity, van der Waals forces, electric field, etc. were considered, they 



would need to be accounted for in (|3.2p . It is obvious from the following ad absurdum 
argument valid for the interface at equilibrium. Indeed, if some forces were omitted from 
(|3.2[) but present in the force balance, (|3.ip would result in the Ap and T % variations along 
the interface. However according to the equilibrium condition, both temperature and 
pressure need to be homogeneous. The T l variation can appear only when the evaporation 
creates the viscous pressure drop due to the fluid flow as shown below in (|4.10|) . 

In principle, the boundary condition (|3.1|) is sufficient to relax the CL singularity and 
simulate the bubble growth on the isothermal heater. However the heat flux q % L (and 
thus the bubble growth rate) would be then overestimated because in reality the heater 
temperature in the vicinity of the CL varies strongly, see sec. 15.51 In a more realistic 
modeling, the heat conduction in the heater needs to be taken into account so that both 
v apour-li q uid and liquid-solid i nterfa ce temperatures allowed to vary. 



Morrisl (|2000h : 



Aiaev et al. 



(|2002l ) assumed the stationary temperature distribution 
inside the heater. However, it is never stationary because the characteristic heat diffusion 
length \/ast (where t is the bubble growth time and a$ is the heater material temperature 
diffusivity) is usually comparable to a characteristic system size. E.g. for the bubble 
growth problem this size is the bubble radius the dynamics of which is controlled by the 
heat supply from the heater. Another reason is that the temperature distribution in the 
heater cannot accommodate instantaneously the changes in geometry caused by the CL 
motion. Below we study a transient heat conduction problem in a semi-infinite heater 
with an arbitrary distribution of the heat flux qs(x,t) at the heater surface, qs will be 
determined later on by matching with the liquid domain solution. 

The energy is supplied to the heater via a homogeneous volume heating (by the electric 
current). In the framework of the present approach any time dependence of the volume 
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heating power can be treated. The heating power per unit volume is chosen in the form 
C/yi to model a heating impulsion at t — 0; C is an arbitrary con s tant. Although by 
other reasons, the same form has been chosen by iNikolavev et al. 



(|200lD . As in that 



article, the reference heat flux can be defined via C as 



q re f = Cy/nask L /(k S y/aL + k Ly /a^), 



(3.3) 



where k means the heat conductivity of the liquid (fcjj or solid (ks)- 

The homogeneous initial temperature Tg(t — 0) = T is chosen as the initial condition. 
The resulting temperatu re of the heater surface obtained with the Green function method 



(jCarslaw fc Jaeger 



19591 ) is then 

T s (x,t)=T + ^CVt--\ 

k s 2iTks J 

The Tn value 



dr 



Qs(x',t ) 

t-T 



■ exp 



IY2 



(x - x') 
Aa s {t - t) 



T = T sat 1 + 



dx'. (3.4) 



(3.5) 



is chosen to be qual to the saturation temperature corrected for the initial pressure jump 
Apo appearing due to the initial interface curvature. The following assumptions are made 
to solve the conjugated problem in the whole domain. 

• The temperature distribution across the liquid film is assumed to be stationary. 
This approximation is valid when the film thickness is smaller than the thermal diffusion 
length \Ja.i,t. The heat flux is then constant across the liquid and 

T s {x)-T l (x) 



Qs(x) = ql(x) = k L 



h(x) 



(3.6) 



• The vapour is assumed to be non conductive. The heat flux to the vapour domain 
is then neglected at the part of the heater surface in contact with vapour. The infinite 
coordinate integration limits in (|3.4| needs to be replaced by the liquid-solid contact area 
Qls — ^Ls(t) because q$ vanishes at the rest of the heating surface. 



Dynamics of the triple contact line on the non-isothermal heater 11 
By using (|3-H 13. 2\ I3.6|) , one obtains the expression 



TT , , x , r (3-7) 

where all the variables are taken at the same time moment. Equalising with (13. 4[) leads 
to the following integral equation for qs(x,t): 



Qs{x,t) 



R' 



h(x, i) 



1 



At 



Qs(x',r ) 

t-T 



■ exp 



A2 



(a; - x 1 ) 
4a s (t - t) 



da:', (3.8) 



27rfc s Jo Jn LS (r) 
which needs to be solved only for x £ FlLs(t)- 

An important consequence of this quite general heat transfer model consists in the 
following: the pressure jump Ap should be finite everywhere, otherwise the temperature 
(|3.7|) will be infinite, which is non-physical. 



4. Interface shape determination in the microregion 

4.1. Relaxing the hydrodynamic CL singularity 

Generally speaking, the interface shape is determined from the solution of the equations 
of hydrodynamics in the liquid and vapour domains with the appropriate boundary con- 
ditions. The viscosity and the inertia of the vapour can be neglected with respect to those 
of the liquid. The equations of hydrodynamics in the vapour then lead to the homoge- 



neous vapour pressure py 



Stephan fc Hammer 



1985 



Analogously t o the thin films treatment (jWavner et al 



Hocking 



1995; 



Anderson fc Davis 



1976 



1995h . the hydrodynamics 



equations in the liquid in the vicinity of the CL can be solved using the lubrication 
approximation. The jump of the normal stress reduces to the pressure jump Ap at the 
interface because the normal viscous stress is negligible in the lubrication approximation. 
For solving the lubrication equations, the boundary condition for the tangential velocity 
v x at the solid surface is necessary to be defined among others. This condition turns out 
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to be extremely important when the CL is allowed to move. It is well known that the 
conventional no-slip condition 

v x = Q (4.1) 

leads to a non-integrable divergence of the stress at the CL so that the force balance 
(see (|4.11l) below) cannot be satisfied. The simplest (and for this reason used by many 
researchers) method of relaxing this singularity consists in using instead of (|4.1ll the 
Navier slip condition 



U™ la 



dv x 
dy 



(4.2) 



Lauga et al 



(120071) . A 



that involves the slip length l s , which is reviewed in detail by 
conventional approach that makes use of the lubrication approximation (see Appendix 
[A"|l shows that the corresponding to (|4.2j) pressure jump is defined by the equation 



d_ 

dx 



h z - + 1, 



dAp 

dx 



PlH 



(4.3) 



wh ere v % = v' l (x) i s the interface velocity. The case q l L — of (|4.3|) has been derived 



by 



Snoeijer et al. 



(J2007|) to describe the contact line motion with no phase change. 



The case q l L = 0, l s = is conventional for the description of the motion of the 



continuous thin films. Another limit 



dl976|); 



Stephan fc Hammerl (jl985l ): 



0,h 



Hockinel (jl995h : 



has been used by 



Anderson fc David (|l995l ). A sim 



Wavner et al 



ila r to (|4.3I) approach (d erived from more general Stokes equation) has been discussed 



by 



Mathieu et al 



(|2002h who also aimed to describe the contact line evaporation. They 
however considered l s — 0; the CL singularity has not been solved. It should be noted 
that q l L and v % are not independent in the case considered her e: the interface moves due to 



evaporation. Similarly to the approach of 
to vary along the interface. 



Aiaev et al 



((2002J) , both q\ and v % are allowed 



Before starting the numerical calculation, an asymptotic approach needs to be per- 
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formed at x — > xcl to determine whether divergencies are encountered. Since q L {x) = 

qsix) is required to be bounded (by the same reasons as Ap), the r. h. s. of (|5.6|) needs to 

be considered constant in this limit because of the finite value of the contact line velocity. 

Close to the CL, 



h ~ 9{x - x C l) 
and the first integration of (|4.3p results in 



3 J ~dx~ ~ ^ ~ CL >' 



(4.4) 



(4.5) 



The subtraction of xcl (which is the integration constant) in the r. h. s. is necessary to 
match the zero value of the 1. h. s. at x = xcl- By using (I4.4[) in (|4.5[) . one arrives to the 
asymptotic expression 

(4.6) 

OX X — XCL 

which results in Ap ~ log \x — xcl\- Being integrable, this divergence does not pose any 
conceptu al problems for the conventional d escription of the moving CL with no mass 



transfer (Hocking 



1995 



Snoeiier et al. 



2007 



For the case of the phase transition at the 



interface, the pressure is however required to be finite at the CL because it determines 
the local interface temperature, see (|3.H3.2|) . The Navier slip condition ()4.2|) is thus 
unsatisfactory. 

The method of relaxing of the CL singularity proposed by IShikhmurzaevI (|1997l) does 



provide a finite pressure at the CL. However we shall not use this method here because 

of its high complexity: it involves additional physical assumptions and requires several 

more equations to be solved in addition to the conventional hydrodynamics. We propose 

instead a modification of the slip boundary condition. The physical significance of the 
| This divergence however does lead to poor convergence of the numerical solution of transient 

problems, see sec. 15.41 
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condition ()4.2|) is transparent. It means simply a linear mobility relation between the 
force exercised on the liquid (i.e. the hydrodynamic stress) and the velocity. The first 
choice would be to assume a nonlinear law (the force proportional to some power of the 
velocity). However, it can be shown (see Appendix |B|) that such a boundary condition 
leads to the logarithmic divergence too. We propose to use instead of (|4.2[) the "second 
derivative" boundary condition 



dv x 



Vx ~ ls dy Pls dy 2 



(4.7) 



where f3 is a constant. The second term involves the normal derivative of the tangen- 
tial viscous stress while the first involves the viscous stress itself. This condition is used 



( Maurer et al. 



2003; 



Hadiiconstantinou 



2003 



Lockerbv et al 



20041 ) to describe high ve- 



locity gas flows. Although the velocities are far from being high in the present case, the 
second derivative term is needed to be accounted for because the stresses are large near 
the CL. The stress normal derivative is negative because the shear stress is a decreasing 
function of the distance from the wall. 

A condition that involves the second velocity derivative is not a natural boundary 
condition for the Navier-Stokes equations. Strictly speaking, the higher order (and for this 
reason much more complicated) Burnett equations ought to be solved instead. However 
there are no technical interdictions to u se the condition (|4 .7I) with the Navier-Stokes 



equations and it is ofte n done in practice Maurer et al 



here too. According to 



Hadjiconstantinovj 



)3) 



2003) . This way will be followed 



20031 ). the value (3 = 0.5 fits well the results 



of numerical simulations; it will be used hereafter. 

The pressure jump corresponding to the second derivative slip condition (|4.7j) is defined 
by the expression (see Appendix [Al for its derivation) 



d_ 

dx 



h ( y + hl s + (ill 



dAp 

dx 



H v 



Vl 
PlH 



(4.8) 
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When performing the asymptotic analysis of this equation, one notices that the term 

(3l 2 in the parentheses of the l.h.s. is prevailing and causes the finiteness of dAp/dx when 

X -> X C L, 



dAp 



dx 



JJ_ 

mi 



v 1 (x C l) 



1l( x cl) 



PlH 



(4.9) 



The pressure jump Ap(xcL) is thus finite and l|4.8[) can be used to describe evaporation 
in the microregion. For the numerical treatment of (|4.8p . it is convenient to enforce the 
condition (|4.9[) by integrating (|4.8p once, which results in 



h 2 

h[— + hl s + I3l\ 



dAp 

dx 



= P 



v l {x') 



PlH 



Ax'. 



(4.10) 



4.2. Interface shape equation 



The pressure jump Ap across the interface can be written as (|Son fc Dhir 



1999) 



Ap = Ka — p r , 



(4.11) 



where K is the interface curvature, a is the surface tension and p r = ^(py 1 — Pl~) is 
the differential vapour recoi l pressure which needs to be taken into account when the 



heater drying is considered (|Nikolavev fc Bevsens 
local mass exchange rate, see (|A 101) . 
In 2D case Eq. (|4.11j) reads 

d 2 h 



1999: 



Nikolavev et al 



20011 ); T) is the 



(4.12) 

where u — dh/dx. The boundary conditions at x = xcl for (I4.12p are given by two 
expressions 



h = 0, 
u = tan ( 



(4.13) 
(4.14) 



The contact angle 9 depends only on the materials of the tree phases at contact as given 
by the classical Young formula. Since 9 and the interface slope in the microregion are 



16 



Vadim S. Nikolayev 

y k 



saturated vapour 




Figure 3. Geometry of the model heat transfer problem of evaporating droplet. 

usually small, (|4. 12114. 14j) reduce to 

d 2 h 



h = and u 



at x = xcl- 



(4.15) 
(4.16) 



Note that the interface shape changes with time because of the time variation of the 
pressure terms. The set of equations (|4.10l) (where qs is used for q l L ) and (|3.8I4.15IC ip 
allows the heat transfer in the microregion to be determined. One more remark concerns 
the CL time evolution xcl = xch(t) that obeys the equation 

dx C L v 1 {x C l) 



dt 



(4.17) 



It follows from (|4.16IC 1IC 3p since the condition Ah /At — holds at the CL. 

One more equ ation is necessary for the problem closure. For example in the bubble 



growth problem (jKern fc Stephan 



2003), the matching of the solutions in the microregion 
(where the thin film approximation is applicable) and the rest of the bubble interface 
(macroregion) is required. 



5. Hydrodynamics and heat transfer during the evaporation of a 
shallow droplet 

The drop evaporation was generally studied for two different evaporation regimes. 
The first is a slow drying of a drop in an atmosphere of another non-condensable gas. 
The evaporation is thus limited by the speed of diffusion of the vapour in the gas. 
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A stationary spatial distribution of the vapour concentration is assumed. It re s ults i n 



Ristenpart et al. 



a well known di s tribut ion!! of q l L along the drop interface, see 



Poulard et al 



(120031) : 



(120071) and references therein. The second drop evaporation regime is 
the opposite limit of the strong evaporation limited by the heat delivery rate. This second 
regime of evaporation of a 2D sessile (posed on a heater, see Fig. [3]) liquid droplet in the 
atmosphere of its saturated vapour will be solved in order to illustrate the application 
of the ideas developed above. We use this example simply because it is easier to treat 
than the bubble growth. The droplet is assumed to be so shallow that the lubrication 
approximatio n can be applied to i ts wh ole surface. A similar problem has been already 



considered by I Anderson &: Davisl (1995) in the approximation of isothermal heater and 
evaporation localised exclusively at the CL. These approximations are dropped here. 
Among other differences one can list our omission of Marangoni flow by the reasons 
already discussed (pure fluid) and the assumption of 9 independent of the CL velocity 



based on more recent results of atomistic simulations of 



Hadi iconst ant inou 



199J). Unlike 



Anderson fc Davisl (|1995h . we do not consider here the wetting hysteresis. Note that it 



can be accounted for in a natural way by int roducing surface heterogeneities modeled 



20051 ). 



with the variation of 9 along the solid surface (jNikolavevi l; 

The particularity of the drop geometry (Fig. [3]) consists in the CL position determi- 
nation. Instead of being determined from the matching of micro and macro regions, it 
is defined by introducing an additional equation, that for the (2D) drop volume V. Its 
time variation is defined by the global heat balance equation 



dV 



qs(x)dx, 



(5.1) 



f The assumption of the stationarity leads to the vapour concentration inversely proportional 
to the distance from the drop. Note that since this function is non-integrable, the heat flux 
expression based on it is not entirely justified. 



18 



Vadim S. Nikolayev 



Variable 



Notation Dimensional quantity 



length 


; Is 


d 


time 


t 


d 2 /a s 


velocity 


v z 


as/d 


pressure jump 


Ap 


a/d 


heat flux 


qs 


q = Cy/nas 


interface resistance 




d/kL 



Table 1. Dimensional quantities used to make the variables dimensionless. 



where the assumption (|3.6p has been applied. 

By taking into account the well known expression 



dV 
~dt 



v\x)dx, (5.2) 



Eq. (|5.ip can be rewritten in the form 

qs(x) 



v\x) 



Ax = Q (5.3) 



PlH 

more convenient to use in the present context. 

The initial shape is chosen separately for each of two problems considered below but 
the contact angle is always equal to 9. The initial drop height defines the length scale d 
of the problem. 

5.1. Reduction to dimensionless form 

The dimensional parameters used as the characteristic scales to reduce the governing 
equations to the dimensionless form are shown in Table [T] where the flux value q used 
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for non-dimensionalising purposes can be related to the reference heat flux (|3.3p 



i = *~t(i + r\l¥)- (5 - 4) 

V k L V as J 

Once made dimensionless, the set of equations (|3.8I4.10I4.15[) reduces to 

r ~ - ■ - Nt 

-q s (x, t)[R l + h(x, t)] + -y-[Ap(z, t) - Ap ] = 

F t XCL (5 ' 5) 

\ / dr / qs(x\r)g(x - x',i - r)dx' , 

V * JO J-XCL 

( ''" ' h% + (3hl 2 ) ^ = f [v\x') - N q q s (x')] dx', (5.6) 



3 ' s J dx 

d 2 h 



-2 



where 



n?2 -Ap + N r qj, (5.7) 



gM = ±. ex p[~) ir,..s, 



1959) 



Ant V 4i 

is the Green function for the transient heat conduction problem (|Carslaw fc Jaeger 
and tilde means the corresponding dimensionless parameter. Eq. (|3.6[) is accounted for 
in (GEU). 

The following dimensionless groups are identified 

r = k s /k L , 
q 2 d 



N r = 



N, 



qd 



9 PiW (5.9) 



N K = 



d PL H' 

qd 
ksT S at ' 
fias 



R l can be seen as an effective dimensionless length scale corresponding to the interfacial 
thermal resistance. iV M represents the dimensionless viscous relaxation time (fid/ a). N r 
measures the strength of the vapour recoil relative to the surface tension. The number 
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N q shows the contribution of the heat diffusion in terms of the latent heat transport. iV e 
shows how far the heating drives the system out of thermal equilibrium. And finally, N a 
characterises the contribution of the surface tension to the variation of the local interface 
temperature. — N a /N e is introduced in (|5.5p for the sake of brevity. 

Eqs. (|5.6H5.7p should be supplied with the boundary conditions corresponding to (|4.16[) 
defined at the CL: 

dh 

h = and — = =p0 at x = ±xcl- (5.10) 
ox 



The dimensionless equation (|5.3[) reads 

1 [v i tf)-N q q s (a/)]da/ = 0. (5.11) 



(5.12) 



The time evolution of the drop surface is defined by the velocity expression 

~i - 

di 

that follows from (IC ID. 



5.2. Symmetry considerations 

In the following, the advantage of the drop symmetry will be taken and only a half of the 
drop < x < xcl will be treated. The integration in (|5.5[ 15 . 1 1[) can then be performed 
over a half of the drop base, 



r ~ ~ ■ ~ Nt 

-q s (Z,t)[R l + h(x,t)] + -±[Ap(x,t) - Ap Q ] = 



TV 



t /■xcl(t) 

dr / qs(x',T)[g(x-x',t-T)+g(x + x',t-T)}dx', 



o Jo 



(5.13) 



(i 



L [v'ix') - N q q s {x')] Ax 1 = 0. (5.14) 



An additional symmetry boundary condition 
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is necessary to be denned in the drop centre x = to replace the boundary conditions 

(|5.10p at x = —xcl- The temporal and spatial behaviour of three variables qg, Ap, h and 

the temporal variation of xcl can now be found numerically. 

Remark: When one uses the differential form (|4.8[) of the pressure equation one needs 

one more boundary condition because of its higher order. The condition 

at x = can be employed. One notices that it follows from the equations (|5.61 15. 14[) . 
The dimensionless deviation of T$ from T sat defined as 

Af s = (T s - T sat )/T sat (5.17) 

takes the following form according to (|3.4[ 13. 5p : 



Af s (x) = Ap N a + 2N e 



7T 



(5.18) 

2N e dr qs(x' ,r)[g(x — x' ,i — t) + g(x + x' ,i — r)]dx' . 



'0 JO 

Once (|5.13j) is solved and qs(%) is known for x <E (0, xcl), (|5.18|) allows Ts(x) to be 
determined over the whole heater surface i € (0, oo). The alternative expression that 
follows from (|3.7p is valid only for x € (0,xcl)- 

Af s = ApN a + N e rqs(R i + h). (5.19) 

5.3. Numerical implementation 

The boundary element method (BEM) can be applied to solve (|5.13|) . Since the problem 
is non- linear it will be solved by iterations. The interval (0,xcl) is divided into m 
"boundary elements" (BE) and the interval (0, t) is divided into F equal subintervals 
At with t ~ FAt. In principle, m can be allowed to change with time b ut in the present 



examp le we will keep it constant. Similarly to the BEM described by 



Nikolavev et al 



( 200lh . the variables are supposed to be constant during each of subintervals and on 
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each BE. The nodes are chosen in the centres of the BEs. The abscissas of the mesh 
points (BE's ends) at the time moment F are denoted xf . The abscissa of the i-th node 
and the BE length are denoted respectively xf F = 0.5(xf_ 1 + xf) and df = xf — x f-x- 
The interval (0,xcl) is remeshed at each time step because of the moving boundary. 
The qs value of at j-th node and during /-th time step is denoted q. . The corresponding 



values of Ap and h are denoted respectively p- and h* . The discrete analog of (|5.13|) , 
that holds for i = 1 ... to is 

F m 



2 Qi 



R l + h. 



F(pre) 



Nt ( F A - x 



FAt 



(5.20) 



/=ii=i 

where the value at the previous iteration is denoted with the upper index pre. The 
symmetrised "coefficients of influence" are defined by the expressions 



G 



fF _ 



f 



fAt 

dx I [g{x- 
(/-l)At 



x, FAt -t)+ g( X V 



, FAt - r)]d7 



X fF (x f j _ 1 + xF)- X fF (x f j +x? F ), 



X fF (x) = 

The coefficient \ can be found analytically: 



/At 



da;' / 5(x',FAt-r)dr. 
C/-i)A* 



(5.21) 



(5.22) 



X fF (x) = f x [x, (F-f + l)At] - / x [z, (F - /)At] , (5.23) 



where 



/xOM) 



40F 



and erfc(-) a nd -Ei(-) are the complement ary error function and the exponential integral 



respectively (jAbramovitz fc Stegun 



19721 ). Note that x (x) = f x (x, At) 



The r emaining 



equat ions (|5.61 15 . 7[) are discretised using the finite volume method 



(FVM) (jPatankar 



19801 ) . The discretised first order equation (|5.6p is written for i 
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l...(m-l) 



F F 
Pi+1 ~ Pi 

+ a i 



3=1 



nF _ -"(^-l) h F _ u 



F-l 



At 



At 



A^f 
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(5.24) 



where 



df 



2h. 



F(pre) 



F{pre) 



/3 + hf {pre) l s +f3i 



and the formula (|C 5[) was employed. The discretised second order Eq. (|5.7|) holds for 
i = 2. . . (m — 1). It reads 



,F 

t+1 

b i+l x i 



„nF _nF 



= df p: 



7V r 



F(pre 



(5.25) 



Accordingly to the order of the equations, three boundary conditions (|5.10l I5.15|) are 
needed: 



h F 2 -hi 

~.nF _ ^nF 
h F - h F 



h F m = -e(x. 

4 { Pi + N r 



F(pre) 



„nF _ rpfiF 



= of < V 



N r 



F(pre) 



(5.26) 
(5.27) 

(5.28) 



The contact line position at time step F is found by solving the equation (|5.14p 

uF-l 



3=1 



nF _ x n{F-l) , F _ , t 







(5.29) 



At 3 At 

with Newton's algorithm. The material parameters for water at 10 MPa on the stainless 
steel (Table [5]) are used for the calculations. The equilibrium contact angle is assumed 
to be 9 = 5 ° . The slip length value is uncertain. According to the detailed analysis of 



Lauga et al 



(2007), its value varies from 1 nm to 1 /im depending on wettability and the 



state of the solid surface. Unless mentioned specifically, l s — 10 nm value is adopted. 



5.4. Isothermal drop retraction 

First, one needs to check if the second derivative boundary condition gives reasonable 
results. We perform here only a brief check. Detailed comparison with existing experimen- 
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Description 


Notation 


Value 


Units 


Saturation temperature 


T sa t 


311 


°C 


Thermal conductivity of liquid 


k L 


0.55 


W/(m K) 


Specific heat of liquid 


CpL 


6.12 


J/(g K) 


Mass density of liquid 


PL 


688.63 


kg/m 3 


Mass density of vapor 


pv 


55.48 


kg/m 3 


Latent heat of vaporization 


H 


1.3 


MJ/kg 


Surface tension 


a 


12.04 


mN/m 


Thermal conductivity of steel 


k s 


15 


W/(m K) 


Specific heat of steel 




0.5 


J/(g K) 


Mass density of steel 




8000 


kg/m 3 



Table 2. Material parameters used in the simulation. 



tal data will be made elsewhere. As a test problem, we use the classical drop retraction 
problem which consists in solving the hydrodynamic part of the problem, eqs. (|5. 24115. 29|) . 
in which the heat fluxes qf are zeroed and the drop volume remains constant throughout 
the evolution. 

Initially, a pancake drop shape is chosen. It is supposed to relax to the final circular 
shape. The drop shape evolution is presented in Fig. |4ji. Fig. [5] shows the evolution of 
the contact line position. The influence of the slip length on the CL dynamics is shown. 
It is evident that l s controls the CL relaxation rate: the smaller l s the slower is the CL 
motion. 

As mentioned in the previous section, both Navier and second derivative slip conditions 
may be applied for the isothermal problem. The influence of the boundary condition on 
the CL dynamics is shown for l s = 10 nm. The behaviour is similar although the Navier 
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(a) (b) 

Figure 4. Dynamics during drop retraction at constant volume (no evaporation), (a) Drop 
shape evolution obtained with the second derivative boundary condition (|4. T[) . The arrows show 
the shape evolution directions. Note the difference of length scales on x and y axes, (b) Spatial 
variation of the pressure jump Ap in the semi-log coordinates for the second derivative (/3 = 0.5) 
and Navier (/3 = 0) slip conditions. One notes the pressure finiteness and divergence in the former 
and latter cases, respectively. The influence of the discretisation (d m i n is the minimal grid size) 
on the results is also shown. 

condition results in a slightly slower contact line motion than the second derivative 
condition. The difference in the pressure behaviour is however much stronger. Fig. [4]b 
shows a typical spatial variation of the pressure jump. It remains negative along the whole 
drop contour and decreases near the contact line which corresponds to the growing by 
absolute value curvature and the apparent contact angl e smaller than 9 wh ich is a usual 



2007I ). Comparing 



feature of the receding CL motion without evaporation (jSnoeiier et 
to the Navier condition case, the second derivative condition results in a smaller variation 
of the pressure in agreement with the asymptotic analysis carried out in sec. I4.ll The 
results of the stability of the numerical algorithm with respect to the spatial meshing 
are also shown. One notices that the grid refinement leads to strong correction of the CL 
pressure for the Navier case, while this correction is negligible for the second derivative 
condition case where the pressure is bounded. In practice this means that the stability 
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60 I 
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f (ms) 

Figure 5. Contact line dynamics in isothermal case. Comparison of the evolutions obtained 
with the Navier and second derivative slip conditions calculated for the same slip length and 
the influence of the slip length on the contact line dynamics within the second derivative slip 
model. 

of the numerical algorithm is much better for the second derivative condition than for 
the Navier condition. 

5.5. Drop evaporation with retraction 

It is well known that there are two regimes of drop evaporation: with and without CL 
retraction. The regime of immobile contact line occurs due to its pinning on the solid 
surface defects and will not be considered in this paper. 

The solution of the full problem taking the evaporation into account results typically 
in the drop surface evolution shown in Fig. [6] The drop is initially has an equilibrium 
shape corresponding to the contact angle 8 = 5°. The initial pressure jump Apo is chosen 
accordingly to the drop initial curvature, Eq. (|4.11[) . Than the bulk heating of the solid 
begins and the drop volume decreases until the drop evaporates completely, see Fig. [7l 

One mentions several features of this evolution. While the drop volume decreases al- 
most linearly, the time evolution of other parameters is more sophisticated. Instead of 
a monotonous decrease that one would expects, the drop height grows about 30% be- 
fore starting to decrease. This occurs because of the impulse heating mode chosen here. 
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° 50 10 20 30 40 50 60 70 

x (|im) t (ms) 



(a) (b) 

Figure 6. Dynamics during the drop evaporation, (a) Drop shape evolution for 
q re f = 0.5W/cm 2 . The arrows show the shape evolution directions. The shape in the CL vicinity 
for t = 21.34 ms (marked approximately with a small rectangle) is zoomed in the insert in the 
upper right corner of the plot where the circles show the calculated interface and the straight line 
indicates the asymptote at the CL. (b) Evolution of the contact line position for two different 
values of g re /. 




5 10 15 20 



f(ms) 

Figure 7. Time evolution of the drop volume and height for q re f = 0.5W/cm 2 . 
Because of the high heat input in the beginning of the evolution, the heat flux attains 
rapidly high values at the CL (Fig. [8b) ■ The difference of heat fluxes at the CL and at 
the top of the drop is especially high in the beginning of the evolution. High heat flux 
means high velocity of the fluid supply causing a high pressure drop at the CL reflected 
in Fig. [8^,. The interface curvature (|4.11|) is defined mainly by the pressure jump since 
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the vapour recoil contribution is small in the range of the heat fluxes considered in the 
simulation. Therefore the curvature in the vicinity of the CL follows almost instanta- 
neously the heating impulsion. Since the curvature defines the rate of spatial change of 
the interface slope, a high curvature means that the apparent (visible) contact angle 
becomes larger than 9. This feature is visible in Fig. [5^ where the actual contact angle is 
always equal to 9. However the apparent contact angle increases right from the beginning 
of the evolution. Note that the vapo ur recoil term also leads to th e apparent contact an- 
gle increase as described in detail bv lNikolavev fc Bevsend (|1999j ); the effect thus should 



become even stronger at higher heat fluxes. The time variation of the volume is relatively 
slow because it follows the averaged over the drop surface heat flux. It is evident that 
the increase of the apparent contact angle at almost fixed volume should result in the 
initial increase of h max (Fig. [7|) . This effect should exist as well in 3D and can thus be ob- 
served experimentally. Although t he apparent contact angle increase during evaporation 



is pre dicted by many approaches (Wav ner et 



1995) and has been observed by 



Poulard et al. 



1976 



2003); 



Hocking 



1995 



Hegseth et al. 



Anderson fc Davis 



(2005) among many 



others, the height "jump" at impulse heating was not studied experimentally. 

Accordingly to this picture, two stages of drop evaporation can be identified in Figs. 
[B]and[7] First, the CL moves fast (Fig. [5b, t < 7 ms for q re f — 0.5W/cm 2 ), the apparent 
contact angle and the drop height grow. On the second stage (t > 7 ms), the apparent 
contact angle stabilises, h max starts to decrease and the CL receding speed becomes 
smaller. 

As a matter of fact, the evaporation and interface velocity terms enter the r.h.s. of 
the pressure equation (|4. lOj) with the opposite signs. On one hand, a positive interface 
velocity leads to the reduction of curvature (like in the drop retraction problem, Fig. 
Hfc). On the other hand, the evaporation causes the curvature increase and it is difficult 
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(a) (b) 

Figure 8. (a) Pressure jump and (b) heat flux across the vapour- liquid interface versus 
distance to the contact line in semi-log and log-log coordinates respectively for 

q re f = 0.5W/cm 2 . 

to predict the pressure behaviour a priori. Fig.[8K shows that the latter tendency wins. 
One can see in Fig. [8^, that the curvature near the CL can become orders of value larger 
than in the centre of the drop. In addition, the pressures in the drop centre and at the 
CL are of the opposite signs so that an inflection point exists. It is not resolved in the 
drop shapes of Fig. [6^i because it is situated at less than 1 fim from the CL according to 
Fig. EH- The positive curvature and the appearance of the large apparent contact angle 
is illustrated in the insert to Fig. [5^. where a close vicinity of the CL is presented. Note 
that the pressure is finite at the CL as it should be according to the asymptotic analysis 
of sec. |4~T1 

An other interesting feature already underlined in many studies (ptephan &: Hammer 



1985 



Anderson fc Davis 



1995 



Nikolayev et al 



200 ll ) is the localisation of the heat and 



mass transfer in a close vicinity of the CL. Fig. [SJd shows that the almost all heat flux 
(and thus evaporation) is concentrated on about 1% of the vapour- liquid interface in a 
close vicinity of the CL. Note that unlike the calculation involving the adsorption film 
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Figure 9. Snapshots of the temperature variation along the heater surface for q re f = 0.5W/cm 2 . 
Note the sharp minimum appearing in the close vicinity of the CL. The times corresponding to 
the curves are given in ms. 



( Wavner et al 



1976: 



Stephan fc Hammer 



1985), the heat flux does not vanish at the CL 



and attains there its maximum. The discontinuity of qs(x) thus occurs at the CL. 

The variation of the temperature (|5.18p along the whole heater surface presented in Fig. 
[9] reflects this discontinuity. At small heat fluxes considered here, the absolute variation 
of the heater surface superheating defined as Tg — T sat is of the order of 0.5 K at this 
heat flux. T$ exhibits a minimum at the CL due to the strong latent heat consumption. 
The minimum is very sharp and corresponds to the locally high heat exchange along the 
heating surface so the temperature gradients are huge in this area and correspond to 
the high heat flux q l L (Fig. [BJa). This explains why the temperature variation along the 
heater needs to be accounted for. The sharp temperatu r e min imum at the CL has been 



o btained both in simu l ation s (see e. g 



Theophanous et al. 



Nikolayev et al 



(|200ll) ) and experiments (see e. 



([20021) ) both for drop evaporation and bubble growth in boiling. 
The results of this section have been obtained using the second derivative boundary 
condition (|4.7p . One needs to mention that the case of Navier slip condition could in 
principle be treated numerically. However it would result in the dependence of the results 
on the spatial discretisation because of the pressure singularity. 
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6. Conclusions 

The adsorption (or prewetting) liquid film that has been assumed to cover continuously 
the heater in prior studies of the liquid meniscus evaporation. On the heated surfaces, 
its thickness is found to be smaller than few nanometres under practically important 
conditions. Since the fluid flow in such films is insignificant, the true triple liquid- vapour- 
heater contact needs to be considered. 

Two major difficulties occur while describing the heat and fluid flow in the vicinity 
of the contact line, namely the thermal and hydrodynamic singularities that need to be 
relaxed. The hydrodynamic singularity is especially important because the Navier slip 
boundary condition used conventionally for its relaxation turns out to be inappropri- 
ate when evaporation/condensation is involved. While the viscous stresses are relaxed 
with the Navier slip condition, the pressure remains infinite and causes a non-physical 
divergence of temperature at the contact line. Another, "second derivative" slip condi- 
tion is proposed. It relaxes both viscous stress and pressure singularities. We show that 
for the conventional dynamic contact line problems with no heat transfer it results in a 
behaviour very similar to that obtained with the Navier slip condition. 

High heat fluxes occurring at the contact line create important temperature gradients 
even in the metal heaters with high heat conductivity. This requires considering the heat 
conduction in the heater. Because of fast change in geometry due to the moving contact 
line, the transient heat conduction in the heater needs to be solved. 

On one hand, the shape of the vapour-liquid interface changes in time because of the 
action of various forces. This shape change induces a non-stationary flow in the bulk of 
the liquid. On the other hand, the interfacial evaporation (condensation) creates another 
flow bringing the liquid to the interface. These both flows are coupled and determine the 
contact line position. 
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An approach describing the contact line motion caused by evaporation (or possibly 

condensation) and integrating all the ingredients mentioned above has been developed 

by using the lubrication approximation. It allows solving of the conjugate problems of 

hydrodynamics and heat transfer (liquid and solid domains) in the "microregion" , a 

vicinity of the contact line where the main part of the heat and mass transfer takes place. 

This approach can be used to describe many practical situations like bubble dynamics 

during the boiling, drop evaporation or condensation etc. 

A problem of evaporation of a 2D drop posed on a heated support in the atmosphere 
of saturated vapour has been considered and analysed under assumption of the ideally 
smooth surface (no pinning). Due to the impulse mode of heating, a brief height increase 
of the drop occurs in spite of the drop volume decrease. This effect is a manifestation of 
the localisation of the evaporation at the contact line and associated with the apparent 
contact angle increase. 

The author is grateful to D. Beysens, P. Colinet, D. Jamet, O. Lebague, B. Mathieu, 
G. Gavrilyuk and P. Stephan for helpful discussions. The partial financial support of 
CNES is acknowledged. 



Appendix A. Lubrication theory for the moving contact line with 
heat and mass transfer 

The lubrication theory has been developed independently by|] 



Petrofi" 


(1883 


) and 


Revnolds 



(| 18861 ) . The case of simultan e ous heat and mass transf er has been treated in several pa- 



pers (see e.g. 



Hockine! (Il995h : 



Anderson fc David (jl995l )V However, the lubrication equa- 



tions were written there in somewhat different form inconvenient for the purposes of the 
present study. For the convenience of the Reader, the employed equations are re-derived 
here. 
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For thin fluid layers, the fluid is supposed to move mainly along x axis (Fig. [2|), i.e. 

v x v y . In addition, the v x variation across the layer is assumed to be much larger than 

along it: dv x /dy ^> dv x /dx. The Stokes equations then reduce to: 

dp L d 2 v x 

(A1) 

^=0. (A2) 
dy 

By taking the y derivative of (|A 1[) and using (|A 2[) one arrives at the equation d 3 v x /dy 3 = 
0, the solution of which is 

v x — a + by + cy 2 , (A3) 

where a, b, c are independent of y. They are to be determined from the boundary condi- 
tions. The first of them reflects the absence of the tangential stresses at the free vapour- 
liquid interface y = h(x) 

<"» 

The volume flux Q flowing through the film at a given position x 

rh 

v x (y)dy (A 5) 

serves as the second equation for three unknowns a, b, c. The third condition is given 
by either (|4.2[) or (|4.7p . The back substitution of the solution into (|A 1|> written at the 
vapour-liquid interface results in the following expressions : 

«~KT +rt -)tr 

and 

Q = -^(Y + h% + m °)l£ (A7) 

written for the boundary conditions l|4.2p and (|4.7|) respectively. 

By using the fluid mass conservation, Q can also be expressed via the component v n 
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of the liquid velocity normal to the vapour-liquid interface: 

rl 
/o 

can be rewritten as 



Q= [ v n (l)dlfa f v n (x)dx. (A! 

JO Jxnr, 



*-£■ < A9 > 

where v n is related to the heat flux at the interface q l L (assumed positive when the heat 
flows from the liquid side, i.e. at evaporation) known from the heat transfer problem via 
the mass conservation law 

V =(v l -v n ) PL =ql/H. (A 10) 

The normal interface velocity v % is considered to be positive if directed inside the liquid 
(as the vector ft in Fig. [2]). It is assumed here that the heat conduction in the vapour is 
negligible with respect to that in the liquid. 

By injecting (|A 10[) and (|A 6p or (|A 7p into (|A 9[) . one arrives finally to the expressions 
(j4~3)) or (|4~T0| obtained with (|4~2)) or (|4~7)) respectively. 



Appendix B. Asymptotic analysis for the nonlinear slip boundary 
condition 

A nonlinear power slip law 

«—(£)' < B1 > 

that generalises (14. 2[> is applied here in the asymptotic limit x — > xql where 7 is some 
power. By repeating the procedure described in Appendix \^ one obtains easily the 
equation 

Q-h=( hdAp Y I h39Ap (B2) 
\fi dx J 3fi dx 
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One can reuse the argument that led to (|4. 5|) and get instead an asymptotic expression 

fhdA P y h 2 dA P 

Let us now consider the terms in the l.h.s. If the second of them was prevailing at 
x — > xcLt the divergence would be non-integrable like in the case of no-slip condition. 
Therefore, the first term must be prevailing. The second term can then be neglected and 
it becomes evident that (|B 31) results in (|4.6p independently of the value of 7. 



Appendix C. Numerical calculation of the velocity of the moving 
interface 

Calculation of the normal velocity of the interface from the time evolution of the 
interface shape is somewhat delicate. The expression 



dh 
"dt 



(1 + u") 



2s-l/2 



dh 
"dt' 



(CI) 



Hocking (|l995l )1 requires 



commonly used in the standard lubrication theory (see e.g. 
that x is maintained constant while calculating the time derivative. For the calculation 
purposes, this derivative needs to be replaced by the finite difference 



dh 



h b {x = xf)-h t -\x = xf) 
At 



(C2) 



where the superscript means the time step number and the subscript enumerates the grid 
points. Eq. (IC 2j) involves the height of the interface h F_1 defined at a grid position xf F 
at the next time step. Since the grid points move when the interface moves, one would 
need to interpolate in order to find h F_1 at x = x™ F because h F ~ 1 is known only at the 
node points . In addition, the contact line motion can make impossible finding of 

h F ~ 1 (x — x™ F ) if Xq^ 1 < xf F < x F L , i.e. if the interface at time F — 1 did not exist at 
x = xf F . Both these difficulties can be avoided if the expression for the full derivative of 
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h(x l (t),t) 

dh dx' 1 dh dh 
dt ~ ~dtd^ + ~di 

is applied. Its finite difference approximation is easy to calculate, 

h F (x = x n L F ) - h F ~\x = xf F ~ X) ) _ hf - h 



dt 



By using this expression, one obtains 



At At 



(C3) 



(C4) 



T nF _ n(F-l) , F _ uF-1 

v\x = xf) - uf" 1 - 1 A 7 (C5) 

and both above mentioned difficulties are avoided. 
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